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We analyse numerically the linear stability of the fully developed flow of a liquid metal in 
a square duct subject to a transverse magnetic field. The walls of the duct perpendicular 
to the magnetic field are perfectly conducting whereas the parallel ones are insulating. In 
a sufficiently strong magnetic field, the fiow consists of two jets at the insulating walls and 
a near-stagnant core. We use a vector stream function formulation and Chebyshev col- 
location method to solve the eigenvalue problem for small-amplitude perturbations. Due 
to the two-fold reflection symmetry of the base flow the disturbances with four different 
parity combinations over the duct cross-section decouple from each other. Magnetic field 
renders the fiow in a square duct linearly unstable at the Hartmann number Ha » 5.7 
with respect to a disturbance whose vorticity component along the magnetic field is even 
across the field and odd along it. For this mode, the minimum of the critical Reynolds 
number Rcc ~ 2018, based on the maximal velocity, is attained at Ha ~ 10. Further 
increase of the magnetic field stabilises this mode with Rcc growing approximately as 
Ha. For Ha > 40, the spanwise parity of the most dangerous disturbance reverses across 
the magnetic field. At Ha w 46 a new pair of most dangerous disturbances appears with 
the parity along the magnetic field being opposite to that of the previous two modes. The 
critical Reynolds number, which is very close for both of these modes, attains a minimum, 
Rec « 1130, at Ha « 70 and increases as Rcc « 91Ha^^^ for Ha :S> 1. The asymptotics of 
the critical wavenumber is kc ~ 0.525i7a^^^ while the critical phase velocity approaches 
0.475 of the maximum jet velocity. 



1. Introduction 

Application of a strong magnetic field to a fiow of an electrically conducting fiuid is 
associated primarily with two effects. First, the magnetic field acts on the mean fiow 



profile often creating infiexion points ( Kakutani (1964)[ ), shear layers (Lehnert (1952) | 



and jets (Hunt (1965)1, thus destabiHsing the otherwise stable fiow. Secondly, strong 
magnetic field tends to damp three-dimensional perturbations making them anisotropic, 
aligned with the magnetic field, and to transform them into quasi-two-dimensional struc- 



tures | |Moffatt (1967)1 [Davidson (1995)| . There are fiows, which combine the effect of 



high electromagnetic damping in some fiow regions, high transverse shear in other re- 
gions, such as jets, and moderate stretching along the magnetic field. Instabilities and 
turbulence may strongly affect the transfer of momentum, heat, and mass in such liquid 
metal fiows which are of major importance for various industrial applications ranging 



from metallurgy and semiconductor crystal growth (Davidson (1999) | to the designs of 



fusion reactors with magnetic confinement (Biihler (2007)). 



Here we will be concerned with linear stability of the fully developed, isothermal, mag- 
netohydrodynamic (MHD) fiow in a constant-area square duct with a pair of perfectly 
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Figure 1. Sketch to the formulation of the problem with the base flow proflle for Ha — 100 (a), 
isolines of base flow {y > 0) and electric current lines {y < 0) for Ha = 6 (a; < 0) and Ha = 100 
{x > 0) shown in the respective quadrants of duct cross-section (b) and the base flow velocity 
proflles at y = for Ha = 0, 6, 50, 100 (c). 



electrically conducting and another pair of perfectly insulating walls in the presence of 
a strong magnetic field. The field is parallel to the insulating walls and perpendicular 
to the conducting ones, and such a flow is known as the Hunt's flow (Hunt (1965)). For 
strong magnetic fields this fiow has a pair of characteristic sidewall jets developing along 
the insulating walls while the velocity in the core of the duct is significantly reduced 
(see figure [T]). These effects are due to the pattern of the electric currents, shown in 
figure [Hb) for y < 0. In the core of the duct the electric currents are induced in the 
direction transverse to the magnetic field and, thus, the resulting electromagnetic force 
nearly balances the pressure gradient driving the fiow. At the insulating side walls, the 
electric current turns almost parallel to the magnetic field and, thus, the electromagnetic 
braking force in these regions is significantly reduced. As a result, the applied pressure 
gradient is balanced there mainly by the viscous shear, and the fiow protrudes through 
the magnetic field in thin jets along the sidewalls. In the limit of a very strong magnetic 
field, the jets, which are of fundamental importance for liquid metal blankets in fu- 
sion reactors (see e.g. [Stieglitz et al. (1996)] [Molokov fc Biihler (1994)[ polokov (1993) 1 
carry almost all of the volume fiux in the so-called parallel layers. Such a velocity pro- 



file is highly unstable, as has been confirmed experimentally by Gelfgat et al. (1971) 



and Platnieks & Freibergs (1972) Linear stability analysis of Hunt's fiow has been at- 



tempted by Fujimura (1989) by assuming two-dimensional mean velocity profile and two- 
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dimensional disturbances, both in the mid-plane of the duct transverse to the field. His 
results are of limited interest owing to the three-dimensional nature of both the mean pro- 
file and the disturbances. A three-dimensional linear stability analysis of a single sidewall 
jet has been carried out by Ting et al. (1991) They consider the flow in a rectangular 



duct with thin conducting walls in the presence of a strong transverse magnetic fleld. 
Although Ting et al. (1991) assume the wall conductance ratio to be small, the magnetic 
field is supposed to be so strong that the relative conductance of both Hartmann and 
parallel layers is even smaller than that of the walls. As a result, the induced electric cur- 
rent passes from the core region directly through the parallel layer into the side wall and 
then to close through normal walls back to the core region. Thus, both the sidewalls and 
Hartmann walls are treated by Ting et al. (1991) as effectively well-conducting bound- 



aries. Direct numerical simulation of this fiow has been undertaken by Miick (2000) for 
the Reynolds number significantly above the Hnear stability threshold for the side layers 



predicted by Ting et al. (1991) 



Here, we present the results of the three-dimensional linear stability analysis of Hunt's 
flow in a square duct, which according to Tatsumi & Yoshimura (1990) is linearly stable 
in the absence of a magnetic field. We show that the instability is far more complex 
than predicted by Ting et al. (1991) and Fujimura (1989) using asymptotic theory and 
two-dimensional approximation, respectively. A magnetic field of moderate strength is 
found to render the flow linearly unstable with respect to two pairs of antisymmetric 
streak-like perturbations of the axial velocity concentrated in the middle part of the 
duct. The most dangerous perturbation is essentially 3D with the component of vorticity 
along the magnetic field being even and odd function across and along the magnetic field, 
respectively. This instability is associated with the appearance of two velocity minima in 
the centre of the duct which at stronger magnetic fields develop into the sidewall jets. As 
the magnetic field strength increases, another essentially 3D instability mode with the 
opposite parity across the magnetic field appears. The critical wavelength of both these 
modes exceeds the width of the duct several times even in relatively strong magnetic 
fields. At the same time, the phase velocity strongly correlates with the maximum jet 
velocity. In a sufficiently strong magnetic field, the critical Reynolds number, based on 
the maximum velocity, increases nearly directly with the magnetic field strength. As 
the sidewall jets develop, two new, much more unstable modes appear with the parity 
along the magnetic field opposite to that of two previous modes. The critical Reynolds 
number, which is almost the same for the last two modes, increases at high magnetic fields 
inversely with the side layer thickness while the critical wavelength reduces directly with 
the thickness. 

The paper is organised as follows. In Section 2 below we formulate the problem. Nu- 
merical method is outlined and verified in §3 and numerical results are discussed in §4. 
Section 5 summarises and concludes the paper. 



2. Problem formulation 

Consider a fiow of an incompressible, viscous, electrically conducting Hquid with density 
/9, kinematic viscosity ly and electrical conductivity a driven by a constant gradient of 
pressure p applied along a straight duct of rectangular cross-section with half-width d 
and half-height h subject to a homogeneous transverse magnetic field B. The walls of the 
duct perpendicular to the magnetic field are perfectly conducting whereas the parallel 
ones are insulating. 
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The velocity distribution of the flow is governed by the Navier-Stokes equation 

1 9 1 

9tv+(v-V)v = — Vp + i^Vv + -f, (2.1) 
P P 

where f = j x B is the electromagnetic body force involving the induced electric current, 
which is governed by the Ohm's law for a moving medium 

j = cr(E + vxB). (2.2) 

The flow is assumed to be sufficiently slow so that the induced magnetic field is negli- 
gible with respect to the imposed one, implying the magnetic Reynolds number Rm = 
HQavad ^ 1, where /ig is the permeability of vacuum and vq is the characteristic velocity 
of the flow. In addition, we assume that the characteristic time of velocity variation is 
much longer than the magnetic diffusion time = ^^ad^ that allows us to use the 
quasi-stationary approximation, according to which E — — V0, where cf) is the elec- 
trostatic potential. The velocity and current satisfy the mass and charge conservation 

V ■ V = V ■ j = 0. Applying the latter to the Ohm's law l|2.2p yields 

V^(t) = 'B-u:, (2.3) 

where u; = V x v is vorticity. At the walls of the duct S, the normal {n) and tangential (r) 
velocity components satisfy the impermeability and no-slip boundary conditions, namely 
Vn\g = and Wt-|s = 0. The conditions for the electric current at insulating and perfectly 
conducting walls are j„|^ = and ir\g = 0, respectively. Boundary conditions for the 
current and velocity applied to Ohm's law result in dn4'\s — and = const for at 
insulating and perfectly conducting walls, respectively. 

We employ the Cartesian coordinates with the origin set at the centre of the duct and 
with X, y and z axes directed along the width, height and length of the duct, respectively, 
as shown in flgurelU with the velocity distribution given by v = {u,v,w). The problem 
admits a purely rectilinear base flow with a single velocity component along the duct 

V = (0, 0, u)(a;, y)) which is shown in flgure [ll[a) for Ha — 100. In the following, all 
variables are non-dimensionalised by using the maximum velocity wq and the half-width 
of the duct d as the velocity and length scales, while the time, pressure, magnetic fleld 
and electrostatic potential are scaled by d'^/i^, pw^, B = |B| and wodB, respectively. Note 
that we use the maximum rather than average velocity as the characteristic scale because 
the stability of this flow is determined by the former as discussed in the following. 

Base flow can more conveniently be described using the z-component of the induced 



magnetic fleld b instead of the electrostatic potential (p (Moreau (1990) I . This temporal 
change of variables does not affect the following linear stability analysis which requires 
the base flow proflle but not the electrostatic potential. Then the governing equations 
for the base flow take the form 

V^w + Hadyb = P, (2.4) 
+ HadyW = 0, (2.5) 



where Ha — dBy/a/{pv) is the Hartmann number and b is scaled by po^Japv^ /d . Note 
that the isolines of b represent electric current Hues which are shown in the bottom part of 
flgurelUb) for Ha = 6 and Ha = 100. The dimensionless constant axial pressure gradient 
P, which drives the flow, is determined from the normalisation condition Wmax ~ 1- 
The velocity satisfles the no-slip boundary condition w = at x = ±1 and y = ±A, 
where A = h/d is the aspect ratio, which is equal to 1 for the square cross-section 
duct considered in this study. The boundary conditions for the induced magnetic fleld 
at insulating and perfectly conducting walls are b = {x = ±1) and dyb = {y — ±A), 
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respectively. The base flow is obtained numerically by the Chebyshev collocation method 
which is described and validated in the next section. 

In order to satisfy the incompressibility constraint V • v = for the flow perturbation, 
we are looking for the velocity distribution as v = V x where tp is a vector stream 
function. The vector stream function as the magnetic vector potential A in electrody- 
namics is defined up to a gradient of an arbitrary function which added to ip does not 
change v. In order to eliminate this ambiguity, we impose an additional constraint on if) 

V-ijj = 0, (2.6) 



which is analogous to the Coulomb gauge for A (Jackson (1998) I . This gauge, similarly 
to the incompressibility constraint for v, leaves only two independent components of ip. 

The pressure gradient is eliminated by applying curl to (|2.ip which yields two dimen- 
sionless equations for ip and uj 

dtuj = V^o; - Reg + Wh, (2.7) 
= V^x() + u}, (2.8) 

where g = V x (v • V)v, and h = V x f are the curls of the dimensionless convective 
inertial and electromagnetic forces, respectively and Re = w^d/v is the Reynolds number. 
In fusion blanket applications Re^ 10^ — 10^ while Ha ^ 10"^ — lO'*. 

The boundary conditions for tp and uj are obtained as follows. The impermeability 
condition applied integrally as v • ds = i/j • dl = to an arbitrary area of the wall 
s encircled by a contour I yields V'tI^ — 0. Using this boundary condition, which implies 
•01^= ntf^nlg , in combination with (|2.6p we obtain dni^nls = 0- In addition, the no-slip 
condition applied integrally v • dl = J^uj ■ ds yields Unls — 0. 

We analyse linear stability of the base flow {V", a>, 4>}{x, y) with respect to inflnitesimal 
disturbances in the form of harmonic waves travelling along the axis of the duct 

u,, 0}(r, t) - u, ^}{x, y) + {^P, ^, $}{x, 2/)e■^*+*^■^ 

where fc is a wavenumber and 7 is, in general, a complex growth rate. Upon substituting 
the solution sought in such a form into l|2.7p , l|2.8p we obtain the governing equations for 
the disturbance amplitudes 

7W = VfeW - Reg + Ha^h, (2.9) 
= V|^ + w, (2.10) 

= Vl4>^LOy, (2.11) 

where Vfe = V + ikez. Because of the solenoidity constraint satisfled by a; similarly to tp, 
we need only the x- and y-components of (|2.9p . namely, — —dxy<j> — dyW, hy — —dyy4> 
and 

gx — k^vw + dyy{vw) + dxy{uw) + i2kdy{wiju), (2-12) 
cjy = —k^uw — dxx{uw) — dxy{vw) — i2kdx{'ww) , (2.13) 

where ii = ik^^{dyytpy - k'^ipy + dxyi'x), v = -ik^^idxxi^x - k'^ipx + dxyipy), and w = 
dx'4'y ~ dyipx- The relevant boundary conditions are 

= i^y ^ dx-li^x = dxi^y - dyljjx = LUx ^ &t X ^ ±1, (2.14) 

^ = V'a; = dyipy = dxiJy - dytpx = cjy = at J/ = ±A. (2.15) 
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Figure 2. Accuracy of the base flow at various Hartmann numbers depending on the number of 
terms in the Fourier series solution (top axis) and the number of collocation points for numerical 
solution (bottom axis). 



3. Numerical method 

We solve the problem posed by (|2.4p - (|2.5p and l|2.9p - l|2.1ip with the boundary condi- 
tions (|2.14p . (|2.15p by a spectral collocation method on a Chebyshev-Lobatto grid using 
even number of points in the x- and y-directions given by 2Nx + 2 and 2Ny + 2, respec- 
tively, where Nx,y = 30 • • ■ 55 is used depending on Ha and Re. The convergence of the 
numerical solution for the base flow was validated against the Fourier series solution of 
Hunt (1965) First, we looked at the relative error in the flow rate for a fixed pressure 



gradient. As seen in figure [21 analytical solution for this quantity converges as 0{N~^) 
and requires « 10'^ terms at Ha — lOP for the relative accuracy of ~ 10^^. The numer- 
ical solution for the flow rate shows a faster-than-algebraic convergence rate developing 
at sufficiently high resolution, which is typical for spectral methods. Second, maximum 
error in velocity scaled with respect to the velocity maximum for — Ny decreases as 
0{N~y). For Ha = 10'^, the resolution of x Ny = 50 x 50 ensures the relative accuracy 
in the base flow velocity of about 10~^. The convergence of the linear stability problem, 
for which the resolution of the base flow is necessary but not sufficient, is tested below. 

Because of the double reflection symmetry of the base flow with respect to x = 
and y = planes, small- amplitude perturbations with different parities in x and y 
decouple from each other. This results in four mutually independent modes which we 
classify as (o, o), (o, e), (e, o), and (e,e) according to whether the x and y symmetry of 
il'x is odd or even, respectively. Our classiflcation of modes specifled in table [T] corre- 
sponds to the symmetries I, II, III, and IV used by Tatsumi & Yoshimura (1990) and 



Uhlmann & Nagata (2006) As a result, the problem is broken up into four independent 
problems of different symmetries deflned in one quadrant of the duct cross-section with 
Nx X Ny internal collocation points. This allows us to reduce the size of the matrix in 
the eigenvalue problem, which is derived below, by a factor of 16. For each symmetry, we 
represent (|2.9p . I|2.10p in the matrix form 



7u;o = AqWo + Aiu;i + fo, 
= BoV^o + '^0, 



(3.1) 
(3.2) 
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I II III IV 

■tpx,uJx,v ■■ (0,0) (o,e) (e,o) (e,e) 

w : (o,e) (0,0) (e,e) (e,o) 

ipz,LJz : (e,o) (e,e) (0,0) (o,e) 

i>y,<^y,u,<j> ■■ (e,e) (e,o) (o,e) (0,0) 

Table 1. The (x, j/)-parities of different variables for symmetries I, II, III and IV; e - even, o - 

odd. 



Nx X Ny 

20 X 20 
24 X 24 
28 X 28 
32 X 32 
36 X 36 



(0.9515252, 
(0.9514767, 
(0.9514760, 
(0.9514760, 
(0.9514760, 



-0.02267611) 
-0.02258387) 
-0.02258432) 
-0.02258432) 
-0.02258432) 



Nx X Ny 

20 X 35 
25 X 40 
30 X 45 
30 X 50 
40 X 60 



(0.23219007, 
(0.23237696, 
(0.23239397, - 
(0.23239343, - 
(0.23239274, - 



-0.3204544 x 10"^) 
-0.5640044 X lO""*) 
-0.28703331 x 10"*) 
-0.21204144 X 10"*) 
-0.21981102 X 10"*) 



Table 2. Convergence of the complex relative phase velocity c = i"f/{Rek) of the least stable 
mode of symmetry I for a model base flow w{x,y) = (1 — x^){A'^ — y^) with A = 1, Re = 10* 
and fc = 1 (le ft) and for the non-magnetic d uct flow with A = 5, Re — 1.04 x 10* and k — 0.91 
considered by Tatsumi & Yoshimura (1990) (right). 



where ^pQ and are the values of {iipx,iJy), {oJx,<^y) at the internal collocation points, 
u)i are unknown values of the tangential component of {ujx^iliy) at + Ny at bound- 
ary points; fo stands for the source term in l|2.9p . Aq, Ai, and Bq matrices represent 
collocation approximation of operator with the explicit boundary conditions (|2.14p , 
l|2.15p eliminated. For the unknown boundary values of Wi, we have an extra boundary 
condition dx'ipy — dynjjx — imposed on (i/'a:, V'y) at + Ny boundary points which is 
represented as 

CoV'o = 0. (3.3) 
To obtain a conventional matrix eigenvalue problem for 7, we need to eliminate loi from 
ISTT]) . f33| . Multiplying both sides of dH]) by CoBg ^ we obtain 

CoB^^Aio;! = -CoBo i(Aou;o + fo), (3.4) 

because ((3?2|) . (|331l imply 

7C0B0 = -jCqiPo = 0. (3.5) 

Now, Wi can be expressed in terms of and fo by solving (|3.4p that substituted back 
into (|3.ip results in 

7a;o = Do(Aoa;o + fo), 
where Do = I — Ai(CoBq ^Ai)"^CoBo"^ and I is the identity matrix. Note that fo is linear 
in both and 0o where the latter can be expressed as — Eq ^a;o,y by solving the 
matrix counterpart of (|2.1ip . Eventually, using l|3.2p . we can write fo ~ ^otpo, that leads 
to 

-fiPo = Bo iDo(AoBo - fo)i>o- (3-6) 
This complex matrix eigenvalue problem is solved by the LAPACK's ZGEEV routine. 

Without the base flow {Re = 0), the leading eigenvalues of l|3.6p are real and negative 
except for Nx + Ny eigenvalues which are zero within machine accuracy. These spurious 
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A Rec X 10"* kc Cc 

5 1.043 0.9085 0.2321 

4 1.819 0.8139 0.2042 

3.5 3.650 0.7075 0.1738 

Table 3. The critical Reynolds number Rec, wavenumber kc and phase velocity Cc obtained 
with 30 X 50 collocation points for the instability mode I in the non-magnetic duct flow at various 
aspect ratios A. 



No: X Ny 

20 X 20 
25 X 25 
30 X 30 
35 X 35 
40 X 40 



(0.7579394, 
(0.7579419, 
(0.7579419, 
(0.7579413, 
(0.7579413, 



{Ha = 10) 
-0.3312792 X 10 
-0.3334821 X 10 
-0.3337441 X 10 
-0.3337346 x 10 
-0.3337034 x 10 

X Ny 

40 X 40 
45 X 45 
50 X 50 
55 X 55 
60 X 60 



N^ X Ny 



25 X 25 
30 X 30 
35 X 35 
40 X 40 
45 X 45 
c {Ha = 10^) 
(0.5053194,0.1453188 x 10 
(0.5054035, 0.1421803 x 10 
(0.5053892,0.1416928 x 10 
(0.5053904, 0.1416891 x 10 
(0.5053902, 0.1417051 x 10 



c {Ha = 10^) 
(0.4907420, -0.8028854 x 10"^ 
(0.4907416, -0.8028572 x 10"^ 
(0.4907415, -0.8028551 x 10"^ 



(0.4907415, 
(0.4907415, 



-0.802855 X 
-0.802855 X 



-2\ 



Table 4. Convergence of the complex relative phase velocity c = i"f/{Rek) of the least stable 
mode for the Hunt's flow in square duct at three different Hartmann numbers: 1) Ha = 10, 



Re = 2000, k = 0.8, mode II; 2) Ha = 10^ 
k = 16, mode I. 



Re = 10^ fc = 5, mode I; 3) Ha = 10^ iZe = 3 x lO'' 



eigenvalues are caused by the way the boundary condition l|3.3p is imposed using l|3.5p 
which can also be satisfied by 7 = 0. These zero eigenvalues can easily be identified and 
discarded. Alternatively, they can be shifted down the spectrum by an arbitrary value 70 
when 7oCot/'o is added to the right-hand side of p.4p . Note that this transformation does 
not affect the true eigenmodes which satisfy the boundary condition l|3.3p . However, our 
approach is not completely free of unstable spurious eigenmodes which may appear at 
sufficiently high Re depending on the collocation approximation of inertial terms (|2.12p , 
l|2.13p . Because the collocation differentiation satisfy the product rule approximately 
rather than exactly (Fornberg (1996)1, the discretisation of inertial terms is affected by 
the form in which they are presented. We find that the number of unstable spurious 
eigenmodes is the least when the inertial terms are approximated in the "conservative" 
form given by (|2.12p . I|2.13p . In contrast to the true eigenmodes, the spurious ones are 
numerical artifacts which depend strongly on the number of collocation points. This 
allows us to identify them easily by recalculating the spectrum with + 1 and Ny + 1 
collocation points and retaining only those eigenvalues whose modulus of the relative 
variation is typically less than s = 10~^ — 10"'*, which is subsequently referred to as the 
relative accuracy threshold. Once a true eigenvalue is identified, it can be tracked further 
by its imaginary part without recalculating the spectrum as the control parameters are 
slowly varied. 

The numerical method has been validated using a model base fiow w{x,y) = (1 — 
x^){A^ — y^) with A = 1 ai Re= 10^ and fc = 1 as well as the non-magnetic duct fiow 



with A = b, Re = 1.04 x 10^ and k = 0.91 considered by Tatsumi fc Yoshimura (1990) 
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Figure 3. Spectrum of the complex relative phase velocities c = i'y/(Rek) for all four mode 
types at Ha = 100, Re — 10^, k — 5 obtained with 50 x 50 collocation points and the relative 
accuracy threshold e = 10""^. 



that resulted in the complex relative phase velocity c — ij/ (Rek) for the least stable mode 
of symmetry I (o, o) which is shown in table [2] for various resolutions. For the model flow, 
the increase of the resolution from 20 x 20 to 28 x 28 collocation points results in the 
fast convergence of the leading eigenvalue with the accuracy raising from two to seven 
figures, respectively, which is comparable to the accuracy of the Galerkin method for 
this test problem used by Uhlmann (2004) Similarly fast convergence is obvious also for 
the non-magnetic duct flow with aspect ratio A = 5. Owing to the large aspect ratio 
^ = 5 as well as the high Reynolds number Re — 1.04 x 10*, which for k ~ 0.91 is 
close to its critical value, at least 30 x 45 collocation points are required to obtain the 
phase velocity with 5 accurate flgures, which again is comparable to the accuracy of the 
Galerkin method tested against the same case by Uhlmann & Nagata (2006) Also the 
instability threshold parameters for the aspect ratios A = 5,4,3.5, which are shown in 



table [3] for 30 x 50 resolution, agree well with Tatsumi & Yoshimura (1990) 



As seen in table HI a comparably fast convergence holds also for the complex phase 
velocity of the least stable modes in the Hunt's flow at Ha — 10, 10^, lO'^. Detailed 
numerical results for these instability modes are presented in the next section. A typical 
spectrum of the complex relative phase velocities c is shown in flgure[3]for Ha — 100 close 
to the instability threshold for the least stable modes of type I and III. The eigenvalues 
have been computed using 50 x 50 collocation points and the relative accuracy threshold 
e= 10-3. 

Subsequently, to verify the numerical accuracy of the obtained results we recalculate 
them with the resolution increased by 5 collocation points in each direction. Only the 
results coinciding by at least four leading flgures are retained. For modes I and III, the 
resolution of 35 x 35 ensures the accuracy of at least 5 digits at Ha = 100, while 50 x 50 
resolution is required at Ha = 3 x lO'^. Modes II and IV require only 30 x 30 resolution 
at Ha « 10, whereas 55 x 55 points are required at Ha ~ 400. 
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Figure 4. Marginal Reynolds number (a) and relative phase velocity (b) versus the 
wavenumber for neutrally stable modes of type II. 



4. Results and discussion 



Here we present the results for the flow in a square duct which according to Tatsumi & Yoshimura (1990) 
is hnearly stable in the non-magnetic case. First, we flnd the base flow numerically and 
normalise it with respect to its maximum velocity which is used here as the velocity scale. 
The flow rate over one quarter of the duct, which is also the average velocity, is found to 
vary for Ha ;:g> 1 as 

Q w 1.23i?a"^/2 + 3.87Ha~\ (4.1) 

where both coefficients are obtained by the best flt of the numerical solution. The main 
contribution to the flow rate is due to the side jets whereas the next-order correction is 
due to the core flow. Although the characteristic velocity of the core flow is only 0{Ha~^) 

1/2 

with respect to that of the side layers, its relative contribution to the flow rate is Ha ' 
times larger because the relative thickness of side jets is 0{Ha~^^^). If the flow rate were 
used for the characteristic velocity, the relative contribution of the core flow in the critical 
Reynolds number would be 0[Ha~^^^). In contrast, when the maximum velocity is used 
for this purpose, the correction is only 0{Ha~^) which becomes negligible at a much 
lower Ha than the previous one. This results in a more deflnite asymptotics appearing 
at numerically attainable values of Ha. Note that for Ha = 10^ the relative contribution 
of the core flow to the flow rate is about 10%. Therefore, we have chosen the maximum 
rather than average velocity as the characteristic scale. 

The neutral stability curves for the instability type II plotted in flgure [4] show the 
marginal Reynolds number, which yields zero growth rate (^[7] = 0) of the most un- 
stable mode for the given wavenumber, and the relative phase velocity c = —uj/{Rek) 
at various Hartmann numbers. Here uj = ^[7] is the frequency of the corresponding 
neutrally stable mode. The minimum of the marginal Reynolds number and the corre- 
sponding wavenumber at which it occurs give, respectively, the critical value Rcc and the 
critical wavenumber kc- This wavenumber along with the corresponding phase velocity 
is plotted in flgure [U against the Hartmann number. It is seen in flgure [Sja) that the 
mode of type II, which is the most unstable up to Ha ~ AQ, flrst appears at Ha « 5.7. 
At this Hartmann number, the velocity proflle of the base flow, which is very close to 
those shown in flgures [Ijb) and (c) for Ha = 6, has a minimum at the centre of the 
duct accompanied by two slight maxima at the each side of it. With the increase of the 
magnetic fleld, these velocity maxima develop into the jets locaHsed at the side walls of 
the duct (see flgure [ifc). There are inflection points in the velocity proflle, which imply a 



Linear stability of Hunt 's flow 



11 



1 1 1 i 1 1 

■ i Ila 


1 


' ',<. 1 ' ' . 

// mode - I 

II ■ 

III 

ly 

' / asymptotic 91//fl 








'.lib /' 
















- ' 1 . . 



Hartttiann numbei". Ha 
0.95 



0.9 
0.85 

0.8 
0.75 

0.7 
0.65 

0.6 
0.55 

0.5 
0.45 



-'.Ha 



(C) 



1 . 1 1 1 1 . 


1 ^^^^^^M^ 


',IIa 




rS^.__ift_ 




(b)„, , 


mode - I ■ 

II ■ 

III ■ 

VI 

^ asymptotic 0.525//fl " " 


lo' 


10" lO' 
Haitmann number, Ha 



■T 

mode - I 
11 

m 

IV 

asymptotic 0.475 



Hartmann number, Ha 



Figure 5. Critical Reynolds number (a), wavenumber (b) and relative phase velocity (c) 

against Hartmann number. 



possibility of an inviscid-type instability, however this criterion is generally restricted to 
one-dimensional inviscid flows (Bayly et al. (1988)). 

Note that in a certain range of the Hartmann number there may be two local minima 
on the neutral stability curve. These are denoted as (Ila) and (lib) in flgurelDJa). The 
flrst minimum, Ila, is below the second one up to « 7 where the critical mode mode 
switches to lib. The corresponding branches of the critical parameters for this mode are 
labelled as Ila and lib in figure [H With the increase of Ha, Rcc first steeply decreases 
down to its minimal value of Rcc « 2018 at Ha ~ 10 and then starts to increase with 
the rate becoming nearly proportional to Ha for Ha > 40. It is important to note that 
the relative phase velocity of the neutrally stable modes, shown in figure [DJb), is nearly 
invariant with wavenumber, and has the order of magnitude 0(1). Moreover, the relative 
phase velocity is seen in figure[5l[c) to stay about 0(1) at large Ha as well. Both of these 
facts imply that the phase velocity of unstable modes is strongly correlated with the 
maximum velocity defined by Re. 

In order to visualise the three-dimensional velocity field of the critical perturbation 
given by 9fi[v(x, ?/)e*'^'=^] we consider the complex amplitude of the velocity perturbation 
V = Vfe X which is associated with the corresponding vector stream function ip. The 
velocity field {u, v) in the {x, j/)-plane can be decomposed into solenoidal Vs and potential 
Vp components which satisfy V • = and V x Vj, = 0, respectively. The solenoidal 
component satisfying the impermeability boundary condition is given by Vs = —e^ x VV'z 
which implies that V'z is the stream function of v^- The incompressibility constraint of 
the whole velocity perturbation, which may be written as V • Vp = ~ikw, in turn, links 
the potential component Vp to the longitudinal velocity perturbation w, which serves as 
a source or a sink for the former. Therefore, the whole velocity perturbation is completely 
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Figure 6. Amplitude distributions of real {y > 0) and imaginary {y < 0) parts of w (a; < 0) and 
lipz {x > 0) of the critical perturbations over one quadrant of duct cross-section for instability 
modes Ila (a) and lib (b) at Ha ~ 7 and Re ~ 10^. 

defined by ipz and w.lna similar way, thie streamlines of solenoidal flow components in the 
(x, z)- and (y, z)-planes are given by i^y and i^x, respectively. Note that the perturbation 
amplitudes are complex quantities whose real and imaginary parts correspond to the 
instantaneous distributions in the (x, j/)-plane shifted in time or in space by a quarter of 
a period. 

Distributions of the most unstable perturbation amplitudes of types Ila and lib are 
plotted in flgure [6] for Ha ~ 7 over different quadrants of the duct cross-section. Both 
perturbations differ mainly by the critical wavenumbers, kc « 0.44 and fcc « 1, respec- 
tively, but have similar amplitude distributions concentrated about the centre of the duct. 
Transversal circulation in the (x, ?/)-plane, which is given by the isolines of ipz, takes place 
about the centre of the duct with Vx and Vy being even functions of x and y, respectively. 
The longitudinal velocity perturbation w, caused by the advection of momentum of the 
base flow by the transversal circulation, is an odd function of both x and y. Both of these 
instability modes are obviously related to the two local velocity maxima which appear 
flrst in the centre of the base flow at Ha « 6. With the increase of Ha these two velocity 
maxima develop into a pair of jets along the insulating side walls (see flgure [TJd, x<0). 

Besides spatial amplitude distributions, the perturbations can be characterised by the 
kinetic energy distribution over the velocity or vorticity /stream function components as 
follows: 

Js Js 

where E is the kinetic energy of perturbation averaged over the wavelength. The inte- 
grals in the expression above are taken over the duct cross-section S, and the asterisk 
denotes the complex conjugate. We flnd that 98% and 91% of kinetic energy for modes 
Ila and lib, respectively, are carried by the longitudinal velocity perturbation w. The 
corresponding component of the vorticity perturbation {uz^ipz), which is associated with 
the circulation in the (x, ?/)-plane, contains only 2% and 7% of kinetic energy, respec- 
tively. The isosurfaces of the critical perturbations of longitudinal velocity and electric 
potential are shown in flgure [7] for one wavelength of mode lib in the right bottom quad- 
rant of duct at Ha = 7 and Re = 10'*. The corresponding perturbation pattern for mode 
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Figure 7. Isosurfaces of longitudinal velocity w (a) and electric potential tp (b) perturbation 
over a wavelength in one quadrant of duct cross-section for instability mode lib at Ha ~ 7 and 
Re « 10''. 
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Figure 8. Marginal Reynolds number (a) and relative phase velocity (b) versus the 
wavenumber for neutrally stable modes of type IV. 



Ila differs mainly by a longer wavelength. As seen in figure [Tlja), the perturbation of w 
represents a pair of elongated, slightly tilted and periodically overlapping streaks located 
close to the centre of the duct. The perturbation of the electric potential, which is the 
largest in the vertical mid-plane of the duct {x — 0), partly reaches the side walls where 
it can be measured experimentally. 

Neutral stability curves for the instability mode of type IV, which appears for Ha > 28 
and differs from the previous one by the opposite x-parity, are plotted in figure [8l Figure 
\5\ shows that Rcc of this mode, which for low values of Ha lies above that of mode II, 
first steeply decreases with Ha by reaching Rcc ~ 9 x 10^ of mode II at Ha w 40. The 
critical Reynolds number for mode IV attains a minimum of Rec ~ 8.8 x 10^ at Ha « 44 
and then starts to increase with Ha remaining below Rcc for mode II up to the largest 
numerically attainable value of Ha « 400. 

Amplitude distributions of the most unstable perturbations of types II and IV are 
plotted in figure [9] at Ha ~ 40 and Re ~ 9 x 10'^. The critical wavenumbers for these 
modes are, respectively, kc ~ 0.75 and kc ~ 0.49. It is seen in figure [9lJ a), that mode II 
has moved from the centre of duct, where it originally appeared at Ha « 5.7, to the side 
wall. The only principal difference between these modes is the opposite x-parity which 
results in a pair of mirror-symmetric longitudinal vortices on each side of the duct with 
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Figure 9. Amplitude distributions of real {y > 0) and imaginary {y < 0) parts of w (a; < 0) and 
{x > 0) of the critical perturbations over one quadrant of duct cross-section for instability 
modes II (a) and IV (b) at Ha « 40 and iZe « 9 x 10*. 

the same or opposite sense of circulation for modes II and IV, respectively. In the first 
case, both vortices are partly connected across the vertical mid-plane of the duct whereas 
they are separated by that plane in the second case. For both modes, the perturbations 
of the longitudinal velocity are localised in the sidewall jets and are very similar to each 
other except for the opposite phases of oscillations across the width of the duct. Modes 
II and IV are also similar from the energetic point of view with 88% and 93% of kinetic 
energy concentrated in the perturbation of the longitudinal velocity. The least amount of 
energy, which is about 1% and 0.4%, respectively, is contained in the x-component of the 
velocity perturbation while the rest carried by the y-component parallel to the magnetic 
field. 

The isosurfaces of the critical perturbations of the longitudinal velocity and of the 
electric potential for mode IV are shown in figure [TO] over one wavelength in the right 
bottom quadrant of duct for Ha = 40 and i?e « 9 x 10'^. In this case, the distributions of iu 
and 4> are even and odd functions of x, respectively. The corresponding pattern for mode 
II at these parameters differs from that of mode IV mainly by the shorter wavelength and 
opposite a:-parity that results in a non-zero perturbation of 4> in the vertical mid-plane 
of the duct (x = 0). 

A pair of additional instability modes of type I and III appears for Ha > 46 and 
47, respectively. These modes differ from the ones of type II and IV by the opposite y- 
parity. The neutral stability curves plotted in figure [TT] look very similar for both of these 
modes. First, for Ha < 54, 60 the neutral stability curves are seen to form closed loops 
which impHes that both modes are unstable only within Hmited ranges of Reynolds and 
wavenumbers. In this range of Ha, there is not only the lower but also the upper critical 
value of Re, by exceeding which all perturbations of the corresponding type become 
linearly stable again. 

These critical values of Re, which are considerably lower than those for the previous two 
modes, are plotted in figure[5]against the Hartmann number along with the corresponding 
wavenumbers and the relative phase velocities. As seen in figure HJa), the upper critical 
Reynolds number steeply increases with Ha becoming very large at Ha « 54, 60 for 
modes I and III, respectively. The lower value of Rec steeply decreases to its minimum 
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Figure 1 1 . Marginal Reynolds number (a) and relative phase velocity (b) versus the 
wavenumber for neutrally stable modes of types I and III. 



Rcc « 1130 attained at Ha w 70. A further increase of Ha results in the growth of the 
critical Reynolds number for both modes approaching the asymptotics Rcc ~ 91Ha^^^ 
for Ha ^ 1. This implies that the critical Reynolds number based on the average velocity 
tends to a constant Rcc ~ 112 while the next-order-correction is about 352Ha~^^^ . In 
contrast to this, the relative next-order-correction for Rcc based on the maximal velocity, 
as discussed at the beginning of this section, is only 0{Ha~^). The critical wavenumber 

1 /2 

for both modes I and II tends to fcc ~ 0.525Ha ' . This means that the critical wavelength 
reduces directly with the characteristic thickness of the parallel layers 0{Ha~^^^). The 
relative phase velocity for both modes is seen to tend asymptotically to a constant Cc ~ 
0.475 which confirms that this instability is indeed associated with the sidewall jets and, 
thus, it is completely determined by the characteristic thickness and by the velocity of 
those jets. Note that the relative phase velocity of two other modes of type II and IV is 
also 0(1) which implies that these instabilities are associated with the sidewall jets, too. 
However, the critical wavenumber for modes II and IV remains 0(1) even for Ha ^ 1. 
This, in turn, implies that both of these instability modes are caused by the velocity 
variation over the height rather than the thickness of the jet. Thus, the height rather 
than thickness of the jet serves as the characteristic length scale for modes II and IV. 
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Figure [12] shows the amplitude distribution of the most unstable perturbation for type 
I only because it is almost identical to that for type III except for the opposite y-parity. 
In this case, however, the y-parity has almost no effect on the amplitude distributions on 
each side of the duct because perturbations are localised in the jets at the side walls and 
practically do not interact with each other. As it is seen, for both modes of type I and 
III, the component of velocity perturbation along the magnetic field is an odd function 
of y whereas the other two velocity components are even functions. Thus, the transversal 
circulation in the (x, ?/)-plane involves a couple of vertically mirror-symmetric vortices 
at each sidewall, while the perturbation of the longitudinal velocity w, which is an even 
function of y, is rather uniform along the magnetic field in the horizontal mid-part of the 
duct. 

From the energetic point of view, it turns out that 70% and 23% of the kinetic energy 
are carried by the z- and x-components of the velocity perturbation, while only 7% 
are carried by the y-component. Thus, 89% of kinetic energy is concentrated in the y- 
component of vorticity /stream function perturbation, which is associated with the z and 
X velocity components. The x- and z-components of vorticity/stream function associated 
with the y-component of velocity contain only 6% and 5% of the energy. Consequently, in 
this case, the perturbation of the fiow is well represented by %py alone whose isosurfaces, 
plotted in figure flSFa) for mode I, show the isolines of solenoidal circulation in the 
horizontal plane. The corresponding isosurfaces of the electric potential perturbation are 
shown in figure [TSlfb) . 

5. Summary and conclusions 

In this study we have analysed numerically the linear stability of the fiow of a liquid 
metal in a square duct subject to a transverse magnetic field. The walls of the duct 
perpendicular and parallel to the magnetic field are perfectly conducting and insulating, 
respectively. We used a novel 3D vector stream function formulation and Chebyshev 
collocation method to solve the eigenvalue problem for small-amplitude perturbations. 
Due to the two-fold refiection symmetry of the base fiow with respect to the a; = and 
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Figure 13. Isosurfaces of ipy (a) and electric potential tp (b) perturbation over a wavelength 
in one quadrant of duct cross-section for instability mode I at Ha ~ 100 and Re ~ 1170. 

y — planes the perturbations with four different parity combinations over the duct 
cross-section decouple from each other. 

The base flow, which without the magnetic fleld is linearly stable in a square duct, 
becomes unstable at the Hartmann number Ha ~ 5.7 as two velocity maxima in the centre 
of the duct appear. This instability mode, which is the most dangerous at low Hartmann 
numbers, involves the vorticity component in the direction of the magnetic fleld which 
is anti-symmetric and, thus, essentially non-uniform along the fleld and symmetric in 
the spanwise direction across the duct. The velocity component in the direction of the 
magnetic fleld for this mode is symmetric along the fleld and anti-symmetric in the 
spanwise direction, respectively. This mode becomes the most unstable at Ha ~ 10 
where its critical Reynolds number based on the maximal velocity attains a minimum of 
Rcc ~ 2018. The increase of the magnetic fleld results in the stabilisation of this mode 
with Rcc growing approximately as Ha. For Ha > 40 another mode with the opposite 
spanwise parity across the duct becomes the most dangerous and remains such up to 
Ha fa 4:6. These two instability modes have most of their kinetic energy concentrated in 
long, streak-like perturbations of the streamwise velocity which flrst appear close to the 
centre of the duct and then move to the sidewall layers as the magnetic fleld increases. 
The critical wavenumber is 0(1) that corresponds to the critical wavelength considerably 
exceeding the width of the duct. It is important to note that the critical phase velocity 
remains 0(1) even in a relatively strong magnetic fleld. This implies that for Ha ^ 1 
both of these instability modes are associated with the sidewall jets. 

At Ha « 46 a pair of two additional instability modes appears with the parity along the 
magnetic fleld being opposite to that of the previous two modes. The critical Reynolds 
number, which is very close for both modes, attains a minimum of Rcc ~ 1130 at Ha « 70 

1 /2 

and increases as Rcc « 91i?a ' for Ha ^ 1. The corresponding critical wavelength is 
kc ~ 0.525Ha^^^ while the critical phase velocity approaches 0.475 of the maximum jet 
velocity. This again suggests these two instability modes, similarly to the flrst two, to be 
associated with the side jets. The main difference between the flrst and second pairs of 
disturbances is in their critical wavenumbers which are 0(1) and 0{Ha^^^), respectively. 
The latter means that the critical wavelength scales directly with the side layer thickness 
0{Ha-^^^), which serves as the characteristic lengthscale for the last two instability 
modes. The critical wavenumber of the flrst two instability modes being 0(1) implies 
that they are associated with the velocity variation over the height of the parallel layer 
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whereas the last two modes are associated with the velocity variation over the thickness 
of this layer. From the energetic point of view, the last two instability modes have most 
of their kinetic energy concentrated in the vortical flow component along the magnetic 
fleld which corresponds to the fluid circulation in the planes transverse to the fleld. 



These last two modes are analogous to the side-layer instability mode found by Ting et al. (1991) 



for the flow in the duct with thin but relatively well-conducting walls, for Ha ^ 1. The 
critical Reynolds number based on the average velocity for the latter is Rcc ~ 313 com- 
pared to our result Rcc « 112 which is rescaled by the average velocity (|4.ip . On the other 



hand, our Rcc is several times higher than the corresponding result of Fujimura (1989) 



for the two-dimensional approximation. Note that this approximation incorrectly pre- 
dicts the base flow in a square duct to remain linearly unstable in the limit of vanishing 
magnetic fleld strength whereas a signiflcant destabiHsation is predicted at the Hartmann 
numbers as small as Ha ~ 10. In addition, note that the instability predicted at Ha w 10 
by our analysis is essentially 3D, as discussed above, and, thus, it is principally differ- 
ent from the 2D one found by Fujimura (1989) Further comparison with the results of 
Ting et al. (1991) shows that our critical wavenumber kc ~ 0.525 scaled by the side-layer 
thickness Ha~^^'^ is close to their asymptotic value kcr = 0.55. At the same time, their 
phase velocity Cc = 0.0947 appears to be signiflcantly lower than ours Cc = 0.423, when 
rescaled with respect to the average velocity. Moreover, the instantaneous streamlines in 
the horizontal mid-plane for the critical perturbation plotted in flgure fTsF a) show dis- 
connected sub- vortices at the sidewall whereas those of Ting et al. (1991) although being 
similarly deformed are fully connected single vortices. These differences may be due to 
the different physical model used by Ting et al. (1991) as discussed in the Introduction. 

In conclusion, note that transiently growing small-amplitude perturbations may appear 
below the linear stability threshold due to the so-called non-normality of the linearised 
operator (Trefethen et al. (1993) I. Transient growth is sought to account for the bypass 
transition to turbulence in the shear flows with a high or none at all Hnear stability 
threshold (Grossmann (2000) I. Such a subcritical transition can hardly be relevant for 
the Hunt's flow in strong magnetic flelds {Ha > 50) because the local critical Reynolds 
number based on the thickness of side layers is already very low ~ 100. However, it may 
still be relevant for weaker magnetic flelds, in which the linear stability threshold is much 
higher or absent at all, when Ha < 5.7. But even in the latter case, linear transient growth 
mechanism might be of Hmited importance because, as argued by Waleffe (1995) '...the 
question of transition is really a question of existence and basin of attraction of nonlinear 
self-sustaining solutions that have little contact with the nonnormal linear problem.^ For 
a non-magnetic square duct flow, such nonlinear self-sustaining solutions in the form of 
flnite amplitude travelling waves have been found recently by Wedin et al. (2009) and 



for the magnetic case by Kinet et al. (2009) 



For a flow in the duct with thin conducting walls, the critical Reynolds numbers 
observed experimentally by Reed & Picologlou (1989) appear considerably higher than 
those predicted by the linear stability theory of Ting et al. (1991) This may be owing to 
the fact that the flow in the experiment was developing with jets accelerating which would 
render them more stable, or that the probes could not reach the thin parallel layer where 
the instabilities occur flrst. However, this may also imply that the side-layer instability, is 
supercritical. Then the delay of the transition to turbulence signiflcantly above the linear 
stability threshold can be accounted for by the distinction between the convective and 
absolute instabilities. The conventional stability analysis presented in this paper yields 
the convective instability threshold at which the flow becomes able to amplify certain 
externally imposed perturbations (Landau & Lifshitz (1987) |. For a small-amplitude su- 
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percritical perturbation in the form of travelling wave to become self-sustained absolute 
instability is necessary (Lifshitz & Pitaevskii (1981) I. 

The authors are indebted to Leverhulme Trust for financial support of this work. 
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